knitr::opts_knit$set(root.dir = '~/Documents/bgs_sims/')
knitr::opts_chunk$set(fig.width=16) 

Libraries and functions

library(tidyverse)
library(ggridges)
library(viridis)
library(cowplot)
orig_theme <- theme_set(theme_cowplot(font_size=20))

Data preparation functions

read_sim<-function(pidata,xidata,tajdata,mu=1.66e-8,neutral=FALSE){
  #read data
  pi<-as.tibble(read.csv(pidata,header=T))
  pi<-mutate(pi,stat=rep("pi",length(pi[,1])))
  xi<-as.tibble(read.csv(xidata,header=T)) 
  xi<-mutate(xi,stat=rep("xi",length(xi[,1])))
  tajD <- as.tibble(read.csv(tajdata,header=T))
  tajD<-mutate(tajD,stat=rep("tajD",length(tajD[,1])))
  #munge into single tibble of means
  data<-rbind(xi,pi,tajD) %>%
    `colnames<-`(c("gen", 1:101, "sim", "stat")) %>%
    gather(window,value,2:102) %>% 
    group_by(gen,stat,window) %>%
    summarize(mean_val=mean(value,na.rm = T)) %>%
    mutate(centimorgen=(as.numeric(window)-51)*20000*mu*100 ,mean_val=mean_val)
    data=mutate(data,mean_val=if_else(centimorgen==0.00,mean_val*1.25,mean_val))
  if(neutral==TRUE){
    data=mutate(data,mean_val=if_else(centimorgen==0.00,mean_val/1.25,mean_val))
  }

  return(data)
}
#neutral=TRUE
#pidata <- "results/tennessen_neutral_pi.csv"
#xidata <-  "results/tennessen_neutral_xi.csv"
#tajdata <- "results/tennessen_neutral_tajD.csv"

Merge neutral and selected values

merge_sel_neutral <- function(sel_sim_data,neut_sim_data){
  merged <- merge(sel_sim_data,neut_sim_data,by=c('gen','stat','window','centimorgen'))
  merged <- mutate(merged,mean_val = if_else(stat=='tajD',mean_val.x-mean_val.y,mean_val.x/mean_val.y))
  names(merged) <- c('gen', 'stat', 'window', 'centimorgen','bgs','neutral','mean_val')
  merged
}

Plot functions

Plot selection landscape

simplot<-function(somedata,mystat,label){
  ggplot(filter(somedata,stat==mystat), aes(x=centimorgen, y=mean_val, group = gen,color=gen)) + 
  scale_color_viridis(name='Generation') +   
  background_grid(major = "xy", minor = "none") +
#  geom_smooth(se=F) + 
    geom_line() +
  ylab(label)  + xlab("Distance (cM)") #+
#  scale_x_continuous(breaks = seq(-1.5,1.5,0.5))
}

simplot_grid <- function(somedata,events,event_names,ratio=TRUE){
  if (ratio == TRUE){
    
    pi_lab <-  expression(bar(pi)/bar(pi)[neut])
    xi_lab <- expression(bar(xi)/bar(xi)[neut])
    d_lab <- expression(bar(D)-bar(D)[neut])
  }
  else{
    pi_lab <- expression(bar(pi))
    xi_lab <- expression(bar(xi))
    d_lab <- expression(bar(D))
  }
legend <- get_legend(simplot(somedata,"pi",expression(bar(pi)/bar(pi)[neut]))+
                       theme(legend.key.height= unit(2, "cm"),
                             legend.title = element_text(size=16),
                             legend.text = element_text(size=11)) +
                       scale_color_viridis(name='Generation',breaks = events, labels=event_names)
                     )
plot_grid(
plot_grid(
simplot(somedata,"pi",pi_lab)+theme(legend.position='none'),
simplot(somedata,"xi",xi_lab)+theme(legend.position="none"),
simplot(somedata,"tajD",d_lab)+theme(legend.position="none"),
ncol = 1,
align = 'hv'),
legend, rel_widths = c(85,15))
}

Function to plot ratios over time and by genetic distance

plot_ratios <- function(data,events){
comb <- data %>%
  filter(window%in%c(51,51-10,51,51+10,51-20,51+20,1,101,51+5,51-5),stat%in%c('pi','xi')) %>%
  mutate(centimorgen=abs(centimorgen)) %>%
  group_by(centimorgen,gen,stat) %>%
  summarise(mean_val=mean(mean_val))
start_val <- comb %>%
  filter(gen==min(comb$gen))

ggplot(comb,aes(gen,mean_val,color=stat)) +
  geom_smooth(method = 'loess',span=0.2)+
#  geom_line(size=1.5) +
  geom_hline(data=data.frame(yint=start_val$mean_val,centimorgen=start_val$centimorgen),aes(yintercept =yint))+
  labs(x=expression(paste('Generation (x ',N[anc],')')),y=expression(x/x[neut]),color='Statistics', 
       caption ='Summary statistics at different distances (cM) from selected region.\nVertical lines mark demographic events. Horizontal lines denote value before demographic changes') + 
#  xlim(-0.1,1.25) +
  facet_grid(.~centimorgen) +
  geom_vline(xintercept = events, size=0.5, alpha=0.2) +
   scale_color_manual(values=c("#999999", "#E69F00"), 
                       breaks=c("pi", "xi"),
                       labels=c(expression(bar(pi)),expression(bar(xi)))) +
  theme(strip.background = element_rect(fill=alpha('lightblue',.2)),
        strip.text.y = element_text(angle = 0,hjust = 0),
        plot.caption=element_text(size=12,hjust = 0)) +
  scale_x_continuous(breaks=c(0,0.3,0.6,0.9,1.2))

}

Function to plot absolute values over time and by genetic distance

plot_over_time <- function(data,events,stats){

  comb <- data %>%
  filter(window%in%c(51,51-10,51,51+10,51-20,51+20,1,101,51+5,51-5),stat==stats) %>%
  mutate(centimorgen=abs(centimorgen)) %>%
  group_by(centimorgen,gen,stat) %>%
  summarise(neutral=mean(neutral),bgs=mean(bgs)) %>%
  gather(selection,value,bgs,neutral) 

start_val <- comb %>%
  filter(gen==min(comb$gen))
stats_label <- c(expression(bar(stats)))
ggplot(comb,aes(gen,value,color=selection)) +
  geom_line(size=1.5) +
  geom_hline(data=data.frame(yint=start_val$value,centimorgen=start_val$centimorgen,sel=start_val$selection),aes(yintercept =yint,color=sel))+
  labs(x=expression(paste('Generation (x ',N[anc],')')),y=stats,color='Selection')+
  facet_grid(.~centimorgen) +
  geom_vline(xintercept = events, size=0.5, alpha=0.2) +
   scale_color_manual(values=alpha(c("#ef8a62", "#67a9cf"),.9), 
                       breaks=c("neutral", "bgs"),
                       labels=c('Neutral','BGS')) +
  scale_x_continuous(breaks = seq(0,2,0.5))+
  theme(strip.background = element_rect(fill=alpha('lightblue',.2)),
        strip.text.y = element_text(angle = 0,hjust = 0),
        plot.caption=element_text(size=12,hjust = 0))
}

Demographies real models

demog_torres <- read.csv('demographies/torres.csv')
demog_torres <- mutate(demog_torres,model='torres')
demog_ten <- read.csv('demographies/tennessen.csv')
demog_ten <- mutate(demog_ten,model='tennessen')
demog_maize <- read.csv('demographies/maize.csv')
demog_maize <- mutate(demog_maize,model='maize')
demog <- rbind(demog_torres,demog_ten,demog_maize)

# Torres events
anc_ext <- 0
OOA_torres <-   0.437017*2
euro_bneck_torres <- OOA_torres + 0.109481*2
torres_events <- c(anc_ext,OOA_torres,euro_bneck_torres)
torres_event_names <- c('Ancestral expansion','OoA bottleneck','Euro bottleneck')
torres_event_sizes <- c(2.10709*18449,2.10709*18449,0.322295*18449)
torres_event_end<- c(1e5,1e5,2e3)

# Tennessen events
OOA_ten <- 0.264196*2
euro_bneck_ten <-OOA_ten+ 0.0763702*2
finel_ext <-euro_bneck_ten+ 0.0502841*2
tennessen_events <- c(anc_ext,OOA_ten,euro_bneck_ten,finel_ext)
tennessen_event_sizes <-  c(1.98*7310,1.98*7310,0.254609*7310,1.339863*7310)
tennessen_event_end<- c(1e5,1e5,1e5,1e3)
tennessen_event_names <- c('Ancestral expansion','OoA bottleneck','Euro bottleneck','Final growth')

# maize events
maize_events <- c(0)
maize_event_names <- c('Domestication bottleneck')

Demography plot

ggplot(demog,aes(x=generation,y=N,color=model))+geom_line(size=2) + 
  scale_y_log10()+ 
  scale_color_discrete(breaks=c("tennessen", "torres","maize"),labels=c('Tennessen et al.','Torres et al.','Bessinger et al.'))+
  # annotate tennessen
  annotate("segment", x = tennessen_events, xend = tennessen_events, 
           y = tennessen_event_end, yend = tennessen_event_sizes, 
           size=1,alpha=0.6,arrow=arrow()) +
  annotate("text", x = tennessen_events,  y = tennessen_event_end,
           label=tennessen_event_names,
           angle=c(90,90,90,0),hjust=0,vjust=c(0.5,0.5,0.5,1)) +
  # annotate torres
  annotate("segment", x = torres_events, xend = torres_events, 
           y = torres_event_end, yend = torres_event_sizes, 
           size=1,alpha=0.6,arrow=arrow()) +
  annotate("text", x = torres_events,  y = torres_event_end,
           label=torres_event_names,
           angle=c(90,90,0),hjust=0,vjust=c(0.5,0.5,1)) +
  # annotate maize
   annotate("segment", x = 0, xend = 0, 
           y = 12278*0.04, yend = 12278, 
           size=1,alpha=0.6,arrow=arrow()) +
  annotate("text", x = 0,  y = 12278*0.04,
           label='Domestication bottleneck',
           angle=0,hjust=0,vjust=1) +
  labs(x=expression(paste('Generation (x ',N[anc],')')),y= "Population size",color='Demography')

ggsave('figures/annotated_demographies.pdf',width = 16)
## Saving 16 x 5 in image

Europeans

Tennessen

Neutral data

tennessen_neutral <-read_sim("results/tennessen_neutral_pi.csv",
                             "results/tennessen_neutral_xi.csv",
                             "results/tennessen_neutral_tajD.csv",
                             neutral=TRUE)

simplot_grid(tennessen_neutral,tennessen_events,tennessen_event_names,ratio = F)
## Scale for 'colour' is already present. Adding another scale for
## 'colour', which will replace the existing scale.

BGS data

tennessen_bgs <-read_sim("results/tennessen_bgs_pi.csv",
                         "results/tennessen_bgs_xi.csv",
                         "results/tennessen_bgs_tajD.csv",
                         neutral=FALSE)

simplot_grid(tennessen_bgs,tennessen_events,tennessen_event_names,ratio = F)
## Scale for 'colour' is already present. Adding another scale for
## 'colour', which will replace the existing scale.

Merge data

tennessen <- merge_sel_neutral(tennessen_bgs,tennessen_neutral)

Plot absolute values

legend <- get_legend(plot_over_time(tennessen,tennessen_events,'pi'))
plot_grid(ncol=2,rel_widths = c(90,10),
plot_grid(
plot_over_time(tennessen,tennessen_events,'pi')+theme(legend.position='none'),
plot_over_time(tennessen,tennessen_events,'xi')+theme(legend.position='none'),
plot_over_time(tennessen,tennessen_events,'tajD')+theme(legend.position='none'),
nrow = 3,align = 'hv',axis = 'lb'),
legend)

ggsave('figures/origVal_distance_summary_tennessen.pdf',width = 16,height = 12)

Plot ratios

simplot_grid(tennessen,tennessen_events,tennessen_event_names,ratio = T)
## Scale for 'colour' is already present. Adding another scale for
## 'colour', which will replace the existing scale.

ggsave('figures/region_summary_tennessen.pdf',width = 12,heigh=9)

plot_ratios(tennessen,tennessen_events)

ggsave('figures/smooth_distance_summary_tennessen.pdf',width = 16,height = 7)

Torres

Neutral data

torres_neutral <-read_sim("results/torres_neutral_pi.csv",
                             "results/torres_neutral_xi.csv",
                             "results/torres_neutral_tajD.csv",
                             neutral=TRUE)

simplot_grid(torres_neutral,torres_events,torres_event_names,ratio = F)
## Scale for 'colour' is already present. Adding another scale for
## 'colour', which will replace the existing scale.

BGS data

torres_bgs <-read_sim("results/torres_bgs_pi.csv",
                         "results/torres_bgs_xi.csv",
                         "results/torres_bgs_tajD.csv",
                         neutral=FALSE)

simplot_grid(torres_bgs,torres_events,torres_event_names,ratio = F)
## Scale for 'colour' is already present. Adding another scale for
## 'colour', which will replace the existing scale.

Merge data

torres <- merge_sel_neutral(torres_bgs,torres_neutral)

Plot absolute values

legend <- get_legend(plot_over_time(torres,torres_events,'pi'))
plot_grid(ncol=2,rel_widths = c(90,10),
plot_grid(
plot_over_time(torres,torres_events,'pi')+theme(legend.position='none'),
plot_over_time(torres,torres_events,'xi')+theme(legend.position='none'),
plot_over_time(torres,torres_events,'tajD')+theme(legend.position='none'),
nrow = 3,align = 'hv',axis = 'lb'),
legend)

ggsave('figures/origVal_distance_summary_torres.pdf',width = 16,height = 12)

Plot ratios

simplot_grid(torres,torres_events,torres_event_names,ratio = T)
## Scale for 'colour' is already present. Adding another scale for
## 'colour', which will replace the existing scale.

ggsave('figures/region_summary_torres.pdf',width = 12,heigh=9)

plot_ratios(torres,torres_events)

ggsave('figures/smooth_distance_summary_torres.pdf',width = 16,height=7)

Generic models

model1

model1_events <- c(0,0.9)
model1_event_names <- c(0,0.9)

Neutral data

model1_neutral <-read_sim("results/generic/model1_neutral_pi.csv",
                             "results/generic/model1_neutral_xi.csv",
                             "results/generic/model1_neutral_tajD.csv",mu=1e-8,
                             neutral=TRUE)
simplot_grid(model1_neutral,model1_events,model1_event_names,ratio = F)
## Scale for 'colour' is already present. Adding another scale for
## 'colour', which will replace the existing scale.

BGS data

model1_bgs <-read_sim("results/generic/model1_bgs_pi.csv",
                         "results/generic/model1_bgs_xi.csv",
                         "results/generic/model1_bgs_tajD.csv",mu=1e-8,
                         neutral=FALSE)

simplot_grid(model1_bgs,model1_events,model1_event_names,ratio = F)
## Scale for 'colour' is already present. Adding another scale for
## 'colour', which will replace the existing scale.

Merge data

model1 <- merge_sel_neutral(model1_bgs,model1_neutral)

Plot absolute values

legend <- get_legend(plot_over_time(model1,model1_events,'pi'))
plot_grid(ncol=2,rel_widths = c(90,10),
plot_grid(
plot_over_time(model1,model1_events,'pi')+theme(legend.position='none'),
plot_over_time(model1,model1_events,'xi')+theme(legend.position='none'),
plot_over_time(model1,model1_events,'tajD')+theme(legend.position='none'),
nrow = 3,align = 'hv',axis = 'lb'),
legend)

ggsave('figures/generic_models/origVal_distance_summary_model1.pdf',width = 16,height = 12)

Plot ratios

simplot_grid(model1,model1_events,model1_event_names,ratio = T)
## Scale for 'colour' is already present. Adding another scale for
## 'colour', which will replace the existing scale.

ggsave('figures/generic_models/region_summary_model1.pdf',width = 12,heigh=9)

plot_ratios(model1,model1_events)

ggsave('figures/generic_models/smooth_distance_summary_model1.pdf',width = 16,height=7)

model2

model2_events <- c(0,0.5)
model2_event_names <- c('Bottleneck','Recovered')

Neutral data

model2_neutral <-read_sim("results/generic/model2_neutral_pi.csv",
                             "results/generic/model2_neutral_xi.csv",
                             "results/generic/model2_neutral_tajD.csv",mu=1e-8,
                             neutral=TRUE)
simplot_grid(model2_neutral,model2_events,model2_event_names,ratio = F)
## Scale for 'colour' is already present. Adding another scale for
## 'colour', which will replace the existing scale.

BGS data

model2_bgs <-read_sim("results/generic/model2_bgs_pi.csv",
                         "results/generic/model2_bgs_xi.csv",
                         "results/generic/model2_bgs_tajD.csv",mu=1e-8,
                         neutral=FALSE)

simplot_grid(model2_bgs,model2_events,model2_event_names,ratio = F)
## Scale for 'colour' is already present. Adding another scale for
## 'colour', which will replace the existing scale.

Merge data

model2 <- merge_sel_neutral(model2_bgs,model2_neutral)

Plot absolute values

legend <- get_legend(plot_over_time(model2,model2_events,'pi'))
plot_grid(ncol=2,rel_widths = c(90,10),
plot_grid(
plot_over_time(model2,model2_events,'pi')+theme(legend.position='none'),
plot_over_time(model2,model2_events,'xi')+theme(legend.position='none'),
plot_over_time(model2,model2_events,'tajD')+theme(legend.position='none'),
nrow = 3,align = 'hv',axis = 'lb'),
legend)

ggsave('figures/generic_models/origVal_distance_summary_model2.pdf',width = 16,height = 12)

Plot ratios

simplot_grid(model2,model2_events,model2_event_names,ratio = T)
## Scale for 'colour' is already present. Adding another scale for
## 'colour', which will replace the existing scale.

ggsave('figures/generic_models/region_summary_model2.pdf',width = 12,heigh=9)

plot_ratios(model2,model2_events)

ggsave('figures/generic_models/smooth_distance_summary_model2.pdf',width = 16,height=7)

model3

model3_events <- c(0,0.6667)
model3_event_names <- c('Bottleneck','Recovered')

Neutral data

model3_neutral <-read_sim("results/generic/model3_neutral_pi.csv",
                             "results/generic/model3_neutral_xi.csv",
                             "results/generic/model3_neutral_tajD.csv",mu=1e-8,
                             neutral=TRUE)
simplot_grid(model3_neutral,model3_events,model3_event_names,ratio = F)
## Scale for 'colour' is already present. Adding another scale for
## 'colour', which will replace the existing scale.

BGS data

model3_bgs <-read_sim("results/generic/model3_bgs_pi.csv",
                         "results/generic/model3_bgs_xi.csv",
                         "results/generic/model3_bgs_tajD.csv",mu=1e-8,
                         neutral=FALSE)

simplot_grid(model3_bgs,model3_events,model3_event_names,ratio = F)
## Scale for 'colour' is already present. Adding another scale for
## 'colour', which will replace the existing scale.

Merge data

model3 <- merge_sel_neutral(model3_bgs,model3_neutral)

Plot absolute values

legend <- get_legend(plot_over_time(model3,model3_events,'pi'))
plot_grid(ncol=2,rel_widths = c(90,10),
plot_grid(
plot_over_time(model3,model3_events,'pi')+theme(legend.position='none'),
plot_over_time(model3,model3_events,'xi')+theme(legend.position='none'),
plot_over_time(model3,model3_events,'tajD')+theme(legend.position='none'),
nrow = 3,align = 'hv',axis = 'lb'),
legend)

ggsave('figures/generic_models/origVal_distance_summary_model3.pdf',width = 16,height = 12)

Plot ratios

simplot_grid(model3,model3_events,model3_event_names,ratio = T)
## Scale for 'colour' is already present. Adding another scale for
## 'colour', which will replace the existing scale.

ggsave('figures/generic_models/region_summary_model3.pdf',width = 12,heigh=9)

plot_ratios(model3,model3_events)

ggsave('figures/generic_models/smooth_distance_summary_model3.pdf',width = 16,height=7)

model4

model4_events <- c(0,0.05,0.5250)
model4_event_names <- c('Bottleneck','Growth','Recovered')

Neutral data

model4_neutral <-read_sim("results/generic/model4_neutral_pi.csv",
                             "results/generic/model4_neutral_xi.csv",
                             "results/generic/model4_neutral_tajD.csv",mu=1e-8,
                             neutral=TRUE)
simplot_grid(model4_neutral,model4_events,model4_event_names,ratio = F)
## Scale for 'colour' is already present. Adding another scale for
## 'colour', which will replace the existing scale.

BGS data

model4_bgs <-read_sim("results/generic/model4_bgs_pi.csv",
                         "results/generic/model4_bgs_xi.csv",
                         "results/generic/model4_bgs_tajD.csv",mu=1e-8,
                         neutral=FALSE)

simplot_grid(model4_bgs,model4_events,model4_event_names,ratio = F)
## Scale for 'colour' is already present. Adding another scale for
## 'colour', which will replace the existing scale.

Merge data

model4 <- merge_sel_neutral(model4_bgs,model4_neutral)

Plot absolute values

legend <- get_legend(plot_over_time(model4,model4_events,'pi'))
plot_grid(ncol=2,rel_widths = c(90,10),
plot_grid(
plot_over_time(model4,model4_events,'pi')+theme(legend.position='none'),
plot_over_time(model4,model4_events,'xi')+theme(legend.position='none'),
plot_over_time(model4,model4_events,'tajD')+theme(legend.position='none'),
nrow = 3,align = 'hv',axis = 'lb'),
legend)

ggsave('figures/generic_models/origVal_distance_summary_model4.pdf',width = 16,height = 12)

Plot ratios

simplot_grid(model4,model4_events,model4_event_names,ratio = T)
## Scale for 'colour' is already present. Adding another scale for
## 'colour', which will replace the existing scale.

ggsave('figures/generic_models/region_summary_model4.pdf',width = 12,heigh=9)

plot_ratios(model4,model4_events)

ggsave('figures/generic_models/smooth_distance_summary_model4.pdf',width = 16,height=7)

model5

model5_events <- c(0,0.05,0.6833)
model5_event_names <- c('Bottleneck','Growth','Recovered')

Neutral data

model5_neutral <-read_sim("results/generic/model5_neutral_pi.csv",
                             "results/generic/model5_neutral_xi.csv",
                             "results/generic/model5_neutral_tajD.csv",mu=1e-8,
                             neutral=TRUE)
simplot_grid(model5_neutral,model5_events,model5_event_names,ratio = F)
## Scale for 'colour' is already present. Adding another scale for
## 'colour', which will replace the existing scale.

BGS data

model5_bgs <-read_sim("results/generic/model5_bgs_pi.csv",
                         "results/generic/model5_bgs_xi.csv",
                         "results/generic/model5_bgs_tajD.csv",mu=1e-8,
                         neutral=FALSE)

simplot_grid(model5_bgs,model5_events,model5_event_names,ratio = F)
## Scale for 'colour' is already present. Adding another scale for
## 'colour', which will replace the existing scale.

Merge data

model5 <- merge_sel_neutral(model5_bgs,model5_neutral)

Plot absolute values

legend <- get_legend(plot_over_time(model5,model5_events,'pi'))
plot_grid(ncol=2,rel_widths = c(90,10),
plot_grid(
plot_over_time(model5,model5_events,'pi')+theme(legend.position='none'),
plot_over_time(model5,model5_events,'xi')+theme(legend.position='none'),
plot_over_time(model5,model5_events,'tajD')+theme(legend.position='none'),
nrow = 3,align = 'hv',axis = 'lb'),
legend)

ggsave('figures/generic_models/origVal_distance_summary_model5.pdf',width = 16,height = 12)

Plot ratios

simplot_grid(model5,model5_events,model5_event_names,ratio = T)
## Scale for 'colour' is already present. Adding another scale for
## 'colour', which will replace the existing scale.

ggsave('figures/generic_models/region_summary_model5.pdf',width = 12,heigh=9)

plot_ratios(model5,model5_events)
## Warning: Removed 2 rows containing non-finite values (stat_smooth).

ggsave('figures/generic_models/smooth_distance_summary_model5.pdf',width = 16,height=7)
## Warning: Removed 2 rows containing non-finite values (stat_smooth).

model6

model6_events <- c(0.9,0.95)
model6_event_names <- c('Bottleneck','Recovered')

Neutral data

model6_neutral <-read_sim("results/generic/model6_neutral_pi.csv",
                             "results/generic/model6_neutral_xi.csv",
                             "results/generic/model6_neutral_tajD.csv",mu=1e-8,
                             neutral=TRUE)
simplot_grid(model6_neutral,model6_events,model6_event_names,ratio = F)
## Scale for 'colour' is already present. Adding another scale for
## 'colour', which will replace the existing scale.

BGS data

model6_bgs <-read_sim("results/generic/model6_bgs_pi.csv",
                         "results/generic/model6_bgs_xi.csv",
                         "results/generic/model6_bgs_tajD.csv",mu=1e-8,
                         neutral=FALSE)

simplot_grid(model6_bgs,model6_events,model6_event_names,ratio = F)
## Scale for 'colour' is already present. Adding another scale for
## 'colour', which will replace the existing scale.

Merge data

model6 <- merge_sel_neutral(model6_bgs,model6_neutral)

Plot absolute values

legend <- get_legend(plot_over_time(model6,model6_events,'pi'))
plot_grid(ncol=2,rel_widths = c(90,10),
plot_grid(
plot_over_time(model6,model6_events,'pi')+theme(legend.position='none'),
plot_over_time(model6,model6_events,'xi')+theme(legend.position='none'),
plot_over_time(model6,model6_events,'tajD')+theme(legend.position='none'),
nrow = 3,align = 'hv',axis = 'lb'),
legend)

ggsave('figures/generic_models/origVal_distance_summary_model6.pdf',width = 16,height = 12)

Plot ratios

simplot_grid(model6,model6_events,model6_event_names,ratio = T)
## Scale for 'colour' is already present. Adding another scale for
## 'colour', which will replace the existing scale.

ggsave('figures/generic_models/region_summary_model6.pdf',width = 12,heigh=9)

plot_ratios(model6,model6_events)

ggsave('figures/generic_models/smooth_distance_summary_model6.pdf',width = 16,height=7)

model7

model7_events <- c(0.9,0.9667)
model7_event_names <- c('Bottleneck','Recovered')

Neutral data

model7_neutral <-read_sim("results/generic/model7_neutral_pi.csv",
                             "results/generic/model7_neutral_xi.csv",
                             "results/generic/model7_neutral_tajD.csv",mu=1e-8,
                             neutral=TRUE)
simplot_grid(model7_neutral,model7_events,model7_event_names,ratio = F)
## Scale for 'colour' is already present. Adding another scale for
## 'colour', which will replace the existing scale.

BGS data

model7_bgs <-read_sim("results/generic/model7_bgs_pi.csv",
                         "results/generic/model7_bgs_xi.csv",
                         "results/generic/model7_bgs_tajD.csv",mu=1e-8,
                         neutral=FALSE)

simplot_grid(model7_bgs,model7_events,model7_event_names,ratio = F)
## Scale for 'colour' is already present. Adding another scale for
## 'colour', which will replace the existing scale.

Merge data

model7 <- merge_sel_neutral(model7_bgs,model7_neutral)

Plot absolute values

legend <- get_legend(plot_over_time(model7,model7_events,'pi'))
plot_grid(ncol=2,rel_widths = c(90,10),
plot_grid(
plot_over_time(model7,model7_events,'pi')+theme(legend.position='none'),
plot_over_time(model7,model7_events,'xi')+theme(legend.position='none'),
plot_over_time(model7,model7_events,'tajD')+theme(legend.position='none'),
nrow = 3,align = 'hv',axis = 'lb'),
legend)

ggsave('figures/generic_models/origVal_distance_summary_model7.pdf',width = 16,height = 12)

Plot ratios

simplot_grid(model7,model7_events,model7_event_names,ratio = T)
## Scale for 'colour' is already present. Adding another scale for
## 'colour', which will replace the existing scale.

ggsave('figures/generic_models/region_summary_model7.pdf',width = 12,heigh=9)

plot_ratios(model7,model7_events)

ggsave('figures/generic_models/smooth_distance_summary_model7.pdf',width = 16,height=7)

model8

model8_events <- c(0.9,0.95,0.9750)
model8_event_names <- c('Bottleneck','Growth','Recovered')

Neutral data

model8_neutral <-read_sim("results/generic/model8_neutral_pi.csv",
                             "results/generic/model8_neutral_xi.csv",
                             "results/generic/model8_neutral_tajD.csv",mu=1e-8,
                             neutral=TRUE)
simplot_grid(model8_neutral,model8_events,model8_event_names,ratio = F)
## Scale for 'colour' is already present. Adding another scale for
## 'colour', which will replace the existing scale.

BGS data

model8_bgs <-read_sim("results/generic/model8_bgs_pi.csv",
                         "results/generic/model8_bgs_xi.csv",
                         "results/generic/model8_bgs_tajD.csv",mu=1e-8,
                         neutral=FALSE)

simplot_grid(model8_bgs,model8_events,model8_event_names,ratio = F)
## Scale for 'colour' is already present. Adding another scale for
## 'colour', which will replace the existing scale.

Merge data

model8 <- merge_sel_neutral(model8_bgs,model8_neutral)

Plot absolute values

legend <- get_legend(plot_over_time(model8,model8_events,'pi'))
plot_grid(ncol=2,rel_widths = c(90,10),
plot_grid(
plot_over_time(model8,model8_events,'pi')+theme(legend.position='none'),
plot_over_time(model8,model8_events,'xi')+theme(legend.position='none'),
plot_over_time(model8,model8_events,'tajD')+theme(legend.position='none'),
nrow = 3,align = 'hv',axis = 'lb'),
legend)

ggsave('figures/generic_models/origVal_distance_summary_model8.pdf',width = 16,height = 12)

Plot ratios

simplot_grid(model8,model8_events,model8_event_names,ratio = T)
## Scale for 'colour' is already present. Adding another scale for
## 'colour', which will replace the existing scale.

ggsave('figures/generic_models/region_summary_model8.pdf',width = 12,heigh=9)

plot_ratios(model8,model8_events)

ggsave('figures/generic_models/smooth_distance_summary_model8.pdf',width = 16,height=7)

model9

model9_events <- c(0.9,0.95,0.9833)
model9_event_names <- c('Bottleneck','Growth','Recovered')

Neutral data

model9_neutral <-read_sim("results/generic/model9_neutral_pi.csv",
                             "results/generic/model9_neutral_xi.csv",
                             "results/generic/model9_neutral_tajD.csv",mu=1e-8,
                             neutral=TRUE)
simplot_grid(model9_neutral,model9_events,model9_event_names,ratio = F)
## Scale for 'colour' is already present. Adding another scale for
## 'colour', which will replace the existing scale.

BGS data

model9_bgs <-read_sim("results/generic/model9_bgs_pi.csv",
                         "results/generic/model9_bgs_xi.csv",
                         "results/generic/model9_bgs_tajD.csv",mu=1e-8,
                         neutral=FALSE)

simplot_grid(model9_bgs,model9_events,model9_event_names,ratio = F)
## Scale for 'colour' is already present. Adding another scale for
## 'colour', which will replace the existing scale.

Merge data

model9 <- merge_sel_neutral(model9_bgs,model9_neutral)

Plot absolute values

legend <- get_legend(plot_over_time(model9,model9_events,'pi'))
plot_grid(ncol=2,rel_widths = c(90,10),
plot_grid(
plot_over_time(model9,model9_events,'pi')+theme(legend.position='none'),
plot_over_time(model9,model9_events,'xi')+theme(legend.position='none'),
plot_over_time(model9,model9_events,'tajD')+theme(legend.position='none'),
nrow = 3,align = 'hv',axis = 'lb'),
legend)

ggsave('figures/generic_models/origVal_distance_summary_model9.pdf',width = 16,height = 12)

Plot ratios

simplot_grid(model9,model9_events,model9_event_names,ratio = T)
## Scale for 'colour' is already present. Adding another scale for
## 'colour', which will replace the existing scale.

ggsave('figures/generic_models/region_summary_model9.pdf',width = 12,heigh=9)

plot_ratios(model9,model9_events)
## Warning: Removed 1 rows containing non-finite values (stat_smooth).

ggsave('figures/generic_models/smooth_distance_summary_model9.pdf',width = 16,height=7)
## Warning: Removed 1 rows containing non-finite values (stat_smooth).